import numpy as np # This is problem 1 def gaussy(A, b, n): s = [0] * n for j in range(0, n): s[j] = abs(max(A[j], key=abs)) l = [i for i in range(n)] for k in range(0, n): ratiomax = 0 maxrow = k for i in range(k, n): ratio = abs(A[l[i]][k] / s[l[i]]) if (ratio > ratiomax): ratiomax = ratio maxrow = i temp = l[maxrow] l[maxrow] = l[k] l[k] = temp for j in range(k + 1, n): multiplier = A[l[j]][k] / A[l[k]][k] for i in range(k + 1, n): A[l[j]][i] = A[l[j]][i] - (multiplier * A[l[k]][i]) b[l[j]] = b[l[j]] - (multiplier * b[l[k]]) #print('This is A', A) #print('This is b', b) x = [0] * n x[n - 1] = b[l[n - 1]] / A[l[n - 1]][n - 1] for j in range(2, n + 1): summ = 0.0 for k in range(n - 1, n - j, -1): summ = summ + A[l[n - j]][k] * x[k] x[n - j] = (b[l[n - j]] - summ) / A[l[n - j]][n - j] print('The solution vector is', x) gaussy([[2, 3, -4, 1], [1, -1, 0, -2], [3, 3, 4, 3], [4, 1, 0, 4]], [-15, 6, 6, -5], 4) print('\n') # This is problem 1a gaussy([[3, -13, 9, 3], [-6, 4, 1, -18], [6, -2, 2, 4], [12, -8, 6, 10]], [-19, -34, 16, 26], 4) print('\n') # This is problem 1b gaussy([[2, 4, -2], [1, 3, 4], [5, 2, 0]], [6, -1, 2], 3) print('\n') # This is problem 1c # n = 6 n = 6 Aa = np.zeros((n, n)) for j in range(0, n): for i in range(0, n): Aa[i][j] = 1 / (i + j + 1) bb = np.zeros((n, 1)) for i in range(0, n): bb[i] = sum(Aa[i]) gaussy(Aa, bb, 6) # n = 15 n = 15 Aa = np.zeros((n, n)) for j in range(0, n): for i in range(0, n): Aa[i][j] = 1 / (i + j + 1) bb = np.zeros((n, 1)) for i in range(0, n): bb[i] = sum(Aa[i]) gaussy(Aa, bb, 15) # This makes sense because when n=15 the condition number is # much larger than when n=6. This means there is much more room # for potential error so it makes sense that the solution vector # doesn't add up to all 1's. # This is problem 2 def seidel(A, b, n, x0, tol): maxiteration = 200 iter = 0.0 abserror = tol + 1 while (abserror > tol and iter < maxiteration): iter = iter + 1 abserror = 0.0 xold = list(x0) for i in range(0, n): summ = 0.0 for j in range(0, i): summ += (A[i][j] * x0[j]) for j in range(i + 1, n): summ += (A[i][j] * x0[j]) x0[i] = (b[i] - summ) / A[i][i] abserror += abs(x0[i] - xold[i]) if (abserror > tol): print('Did not converge') else: print('Converges to', x0) seidel([[3, 1], [2, 6]], [5, 9], 2, [1, 1], 10**(-2)) seidel([[9, -3], [-2, 8]], [6, -4], 2, [0, 0], 10**(-2)) seidel([[1, -3], [-2, 8]], [6, -4], 2, [0, 0], 10**(-2)) # There is a relatively large change in the solutions when you # change just one number. This could be because there's a lot of # room for potential error. This could be due to a possible # large condition number.